function gmm = gmm_fitem(X, num_components)    
    N = size(X, 1);
    dim = size(X, 2);
    
    mu = cell(num_components, 1);
    sigma = cell(num_components, 1);
    mixing_pi = ones(1, num_components);
    clustering = k_means(X, num_components);
    for k = 1 : num_components
        mu{k} = mean(X(clustering == k, :))';
        sigma{k} = zeros(dim, dim);
        for n = 1 : N
            if clustering(n) == k
                x = X(n,:)';
                sigma{k} = sigma{k} + (x - mu{k}) * (x - mu{k})';
            end
        end
        sigma{k} = 1e-5 * eye(dim) + sigma{k} / sum(clustering == k);
        mixing_pi(k) = sum(clustering == k) / N;
    end
        
    gamma = zeros(N, num_components);
    
    last_logp = 1000;
    logp = 1000;
    
    logdets = zeros(1, num_components);
    sigmainvs = cell(num_components, 1);
    for k = 1 : num_components
        logdets(k) = logdet(sigma{k}, 'chol');
        sigmainvs{k} = inv(sigma{k});
    end
    
    fprintf('Running EM for at most 100 iterations ...\n');
    
    for runs = 1 : 100
        
        % E step
        % Calculate loglikelihoods
        loglikelihood = bsxfun(@plus, zeros(N, num_components), log(mixing_pi));
        loglikelihood = loglikelihood - dim/2 * log(2 * pi);
        loglikelihood = bsxfun(@plus, loglikelihood, -1/2 * logdets);
        for k = 1 : num_components
            V = bsxfun(@minus, X, mu{k}');
            loglikelihood(:,k) = loglikelihood(:,k) - 0.5 * dot(V, V * sigmainvs{k}, 2);
        end
        
        % Check if converged
        last_logp = logp;
        logp = sum(logsumexp(loglikelihood, 2), 1);
        
        fprintf('.');
        if (mod(runs, 50) == 0)
            fprintf('\n');
        end
        
        if abs(last_logp - logp) < 0.001
            break
        end
        
        % Calculate gammas
        loglikelihood = bsxfun(@minus, loglikelihood, logsumexp(loglikelihood, 2));
        gamma = bsxfun(@rdivide, exp(loglikelihood), sum(exp(loglikelihood), 2));
        
        % Previous version, not vectorized
        %for n = 1 : N
        %    loglikelihood = log(mixing_pi);
        %    for k = 1 : num_components
        %        loglikelihood(k) = loglikelihood(k) - dim/2 * log(2 * pi) -1/2 * logdets(k);
        %        x = X(n,:)';
        %        loglikelihood(k) = loglikelihood(k) - 0.5 * (x - mu{k})' * sigmainvs{k} * (x - mu{k});
        %    end
        %    loglikelihood = loglikelihood - logsumexp(loglikelihood);
        %    gamma(n,:) = exp(loglikelihood) / sum(exp(loglikelihood));
            
            %for k = 1 : num_components
            %    tmpp = mixing_pi(k) * mvnpdf(X(n,:), mu{k}', sigma{k});
            %    gamma(n,k) = tmpp;
            %    ksum = ksum + tmpp;
            %end
            %gamma(n,:) = gamma(n,:) / ksum;
        %end
        
        % M step - update parameters
        classN = sum(gamma, 1);
        
        for k = 1 : num_components
            mu{k} = 1 / classN(k) * (gamma(:, k)' * X)';
            
            sigma{k} = zeros(dim, dim);
            V = bsxfun(@minus, X, mu{k}');
            sigma{k} = V' * bsxfun(@times, V, gamma(:,k));
            
            % Previous version, not vectorized
            %for n = 1 : N
            %    x = X(n,:)';
            %    sigma{k} = sigma{k} + gamma(n, k) * (x - mu{k}) * (x - mu{k})';
            %end
            sigma{k} = 1e-5 * eye(dim) + sigma{k} / classN(k);
            
            mixing_pi(k) = classN(k) / N;
        end
        
        
        logdets = zeros(1, num_components);
        sigmainvs = cell(num_components, 1);
        for k = 1 : num_components
            logdets(k) = logdet(sigma{k}, 'chol');
            sigmainvs{k} = inv(sigma{k});
        end
        
        %last_logp = logp;
        %logp = 0;
        %for n = 1 : N
        %    loglikelihood = log(mixing_pi);
        %    for k = 1 : num_components
        %        loglikelihood(k) = loglikelihood(k) - dim/2 * log(2 * pi) -1/2 * logdets(k);
        %        x = X(n,:)';
        %        loglikelihood(k) = loglikelihood(k) - 0.5 * (x - mu{k})' * sigmainvs{k} * (x - mu{k});
        %    end
        %    logp = logp + logsumexp(loglikelihood);
            %tmp = 0;
            %for k = 1 : num_components
            %tmp = tmp + mixing_pi(k) * mvnpdf(X(n,:), mu{k}', sigma{k});
            %end
            %logp = logp + log(tmp);
        %end
        
        %logp
        
        %if abs(last_logp - logp) < 0.001
        %    break
        %end
    end
    
    fprintf('\n');
    
    gmm.mixing_pi = mixing_pi;
    gmm.mu = mu;
    gmm.sigma = sigma;
end